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nonhomogeneous  bars  subjected  to  Saint-Venant  torsion  as  a 
variational  problem,  a  necessary  optimality  condition  is 
derived  and  thus  an  iterative  procedure  is  established  based 
on  the  derived  condition.  This  simple  iterative  procedure 
combined  with  the  Finite  Element  Method  (Hybrid  stress 
approach)  is  applicable  to  the  bars  with  different 
cross-sections  other  than  a  square,  although  only  the  square 
cross-section  is  considered  in  this  study. 

From  the  computation  results  performed  by  the  AMDAHL 
470V/6  in  the  Computing  Service,  University  of  Alberta,  it 
is  found  that  the  optimal  solution  in  this  case  where  the 
shear  modulus  function  varies  jump-like  behaves  similarly 
with  that  in  which  the  material  constant  varies 
continuously.  It  is  further  observed  that  the  improvement 
in  the  torsional  rigidity  by  optimization  is  more 
significant  for  a  square  bar  than  for  a  circular  bar. 

Hoping  that  you  find  this  thesis  satisfactory, 


Yours  truly, 


ABSTRACT 


The  growing  interest  over  the  use  of  composite  material 
systems  makes  it  important  to  study  the  problem  of  optimal 
nonhomogeneity  of  a  1 i near ly-e 1 ast i c  bar  in  torsion.  This 
is  a  basic  example  of  the  problem  of  finding  the  most 
efficient  use  of  material  in  such  systems. 

The  optimization  problem  is  formulated  as  a  variational 
problem  in  order  to  derive  the  necessary  conditions  for 
optimality.  In  this  formulation,  the  shear  modulus  function 
whose  values  take  either  of  two  different  shear  moduli  G 
or  G  ,  is  accepted  as  a  control  function  to  be  optimized, 
while  the  proportion  of  rei nforcement  p  is  specified  as  a 
constraint . 

For  the  given  material  constants  G^  ,  G ^  and  a  given 
proportion  of  reinforcement  p,  the  iteration  procedure  is 
carried  out  by  the  computer  program  to  achieve  the  maximum 
(or  minimum)  torsional  rigidity.  The  computer  program  uses 
the  Constant  Strain  Triangle  Finite  Element  Method 
formulated  on  the  basis  of  the  hybrid  stress  approach  for 
calculating  the  stress  solution  of  torsion. 

A  square,  prismatic  bar  of  unit  length  is  taken  for  a 
test  analysis  to  demonstrate  the  correctness  and  accuracy  of 
the  computer  program.  The  test  results  show  that  the 
optimal  solution  is  unique  and  convergent  in  the  practical 
sense . 

Numerical  analysis  is  then  carried  out  for  a  square  bar 
with  different  shear  modulus  ratios  and  with  various 
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proportions  of  reinforcement.  The  results  are  presented  and 
compared  with  the  Known  solutions  for  a  bar  with  circular 
cross-section.  These  analytical  results  show  that  a  square 
bar  can  be  more  efficiently  optimized  than  a  circular  bar. 

It  is  believed  that  the  present  study  may  find 
applications  in  other  fields  such  as  structures,  heat 
transfer  and  seepages. 
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NOMENCLATURE 


A 

area  of  a  cross-section 

A\  'A2 

areas  of  zones  of  material  1  and  2 

a 

radius  of  a  circular  cross-section 

ai 'aj 'ak 

aj=X|<-Xjin  cyclic  order  of  i  ,j  and  k 

b 

radius  of  inner  circle 

C 

boundary  of  a  cross-section 

eff. 

efficiency  of  optimization 

GrG(xry) 

shear  modulus  function 

G1  f  G2 

shear  moduli  of  two  different  materials 

J 

torsional  rigidity 

J-)  r^2  rL^Op 

torsional  rigidities  of  G] -homogeneous ,  Gj- 

homogeneous  and  optimal  nonhomogeneous  bars 

7 

a  side  of  a  square  cross-section 

M 

twisting  moment 

n 

number  of  nodal  points 

P 

proportion  of  reinforcement 

R 

cross-section  of  a  bar 

U 

strain  energy 

u,u(x,y) 

specific  compliance,  which  is  inverse  of  shear 

modulus  function 

U1  fU2 

specific  compliances  of  two  materials 

f  Vy  f  V2 

di spa lcements  in  x,y  and  z  direction, 

respectively 

l/I 

potential  of  the  surface  tractions 

XrYfZ 

rectangular  coordinates 

V 

admissible  functions 
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6 

T 

n 

<t>,<p(xfy) 

^Mxfy) 


twist  angle  per  unit  length 
shear  stress 
complementary  energy 
Prandtl's  stress  function 
warping  function 


Sub sen ipts 

irjfk  three  nodes  of  an  element  numbered  counter  - 

c lockwi se 

1  f2  two  different  materials 

r,r+1  iteration  steps 
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1.  INTRODUCTION 


The  torsion  problem  for  a  prismatic  bar  was  correctly 
formulated  by  Sai nt - Venant [ 20 ] , 1  who  showed  how  the  torsion 
problem  fitted  into  the  general  theory  of  elasticity. 
According  to  Sai nt -Venant ,  the  behavior  of  a  prismatic  bar 
can  be  approximated  by  considering  only  the  behavior  of  a 
typical  cross-sect  ion  of  the  bar,  and  the  twisting  of  the 
bar  is  composed  of  a  rotation  of  this  cross-section  which  is 
proportional  to  its  distance  from  the  origin  and  of  a 
warping  of  the  cross-section  that  is  identical  for  all 
cross-sections . 

The  torsion  function  (warping  function)  which  is  a 
function  of  the  planar  coordinates  only,  is  Known  to  satisfy 
the  two  dimensional  Laplace  equation  thus  the  torsion 
problem  reduces  to  a  Neumann- type  problem.  Alternatively, 
it  is  found  that  the  torsion  problem  can  be  reduced  to  a 
Dir ichlet-type  problem  by  introducing  Prandtl's  stress 
function.  The  solutions  for  torsion  problems  in  terms  of 
these  two  functions  allow  the  complete  analysis  in  the  case 
of  a  homogeneous,  elastic  bar.  For  all  but  certain 
cross-sections,  an  analytical  solution  may  involve  the  use 
of  Fourier  series  or  complex  variables. 

However,  even  these  cumbersome  analytical  solutions  are 
not  generally  available  for  the  cases  of  nonhomogeneous , 
irregular  shape  of  cross-sections  or  bars  in  elastic-plastic 
deformation,  due  to  the  complexity  of  the  analytical 

1  Numbers  in  square  brackets  refer  to  references. 
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procedure . 

These  more  complicated  torsion  problems  have  been  the 
object  of  many  investigations  using  various  techniques  such 
as  a  graphical  method  for  plastic  mater i a  1 s [ 1 4 , 1 9 ]  and 
finite  di f f erence [ 9 ]  or  direct  numerical  technique[4]  for 
the  elastic  case. 

The  Finite  Element  Method  was  employed  by  Herrman[7]  in 
solving  torsion  problems  for  irregular  shapes.  His  work  is 
based  on  the  displacement  formulation,  in  which  the 
potential  energy  is  expressed  in  terms  of  the  assumed 
warping  function,  and  its  stationary  value  leads  to 
approximate  equilibrium  equations.  This  warping  function 
must  be  chosen  such  that  it  satisfies  the  compatibility 
relations . 

In  a  series  of  achievements  by  Moan[15],  Herakovich 
[5,6]  and  Hodge[5,8],  the  finite  element  stress  formulation 
was  used.  Here,  the  stress  function  which  satisfies 
equilibrium  is  substituted  into  the  complementary  energy, 
and  its  stationary  value  yields  approximate  compatibility 
equations.  The  choice  of  "the  best"  element  is  discussed  in 
the  work  of  Moan,  while  quadratic  programming  is  used  by 
Herakovich  and  Hodge  in  solving  elastic-plastic  torsion  of 
hollow  and  nonhomogeneous  bars. 

For  the  assumed  unknown  functions  in  either  the 
displacement  or  stress  formulation,  the  equilibrium  or  the 
compatibility  relations  are  satisfied  only  in  an  approximate 
sense.  For  example,  in  the  case  of  the  displacement 
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formulation,  the  stress  boundary  condition  that  the  normal 
component  of  the  stress  along  the  boundary  of  the 
cross-section  vanishes  is  not  necessarily  satisfied;  while 
in  the  case  of  stress  formulation,  the  warping  is  not 
continuous  over  the  cross-section  especially  in  the  case  of 
a  nonhomogeneous  bar.  In  fact,  the  nonhomogeneous  bar 
studied  in  the  work  of  Herakovi ch [ 6 ]  is  not  truly 
nonhomogeneous  -  the  two  materials  have  different  yield 
limits,  but  the  elastic  shear  modulus  is  constant  throughout 
the  cross-section  of  the  bar,  that  is,  the  bar  is 
nonhomogeneous  only  for  plastic  deformation. 

In  order  to  overcome  some  of  the  deficiencies  of  both 
the  displacement  and  stress  formulation,  the  hybrid  approach 
was  first  introduced  by  P i an [183  and  developed  further  by 
Yamada  et  a  1 . [22] .  A  number  of  variations  are  available 
with  the  hybrid  approach,  however  the  torsion  solutions  for 
nonhomogeneous  bars  can  be  improved  by  assuming  a  stress 
distribution  inside  the  element  and  the  displacements  on  the 
boundaries  of  the  element  so  that  all  of  the  equilibrium, 
compatibility  and  stress  boundary  conditions  can  be 
satisfied. 

In  an  alternate  mixed  approach  which  was  recently 
presented  by  Noor  and  Anderson [ 1 7 ] ,  both  the  displacement 
and  stress  functions  are  assumed  inside  the  element  and  on 
the  boundaries  of  the  element  as  well.  The  solutions 
obtained  by  this  approach  are  much  improved,  but  the  number 
of  equations  for  the  same  number  of  elements  is  larger  than 
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that  which  occurs  in  the  displacement,  stress  or  hybrid 
approaches.  As  a  result  of  these  considerations,  the 
solution  in  this  work  uses  the  hybrid  stress  approach. 

Due  to  the  recent  interest  in  the  development  of 
efficient  composite  material  systems,  the  torsion  problem 
for  nonhomogeneous  bars  has  received  renewed  consideration 
as  an  optimal  design  problem.  The  mathematical  formulation 
of  the  problem  of  finding  the  optimal  control  function 
(shear  modulus  function  in  this  case)  when  the  state 
equation  is  expressed  in  a  partial  differential  form  is  due 
to  Cea  and  Ma 1 anowski [ 2 ] ,  In  their  work,  it  is  shown  that 
the  solution  of  this  type  of  problem  exists  and  is  unique 
and  that  a  convergent  iterative  algorithm  can  be  constructed 
to  find  the  optimal  solution. 

Klosowicz  and  L ur i e [ 10]  developed  the  idea  of  Cea  and 
Malanowski  further  by  employing  a  Lagrange  multiplier  and 
thus  obtained  a  similar  algorithm  to  that  described  in  [2]. 
Klosowicz[9]  then  carried  out  the  iteration  by  using  the 
finite  diffference  technique  for  solving  the  state  equation 
at  every  step  of  the  iteration,  and  finally  obtained  the 
solutions  for  square  and  elliptic  bar  cross-sections.  It 
was  revealed  by  his  study  that  7  the  rigid  reinforcement 
should  be  located  near  the  center  of  the  lateral  surfaces  of 
the  bar  when  the  weak  material  prevails  while  with 
increasing  amount  of  stiffer  material  the  configuration  of 
the  rei nforcement  must  change  to  form  a  belt  around  the 
center  of  the  cross-section7  . 


5 


For  this  same  purpose,  Banichuk[1]  used  a  perturbation 
technique  to  find  the  optimizing  shape  of  cross-section 
instead  of  optimizing  the  shear  modulus  distribution  within 
the  given  cross-section.  Among  all  the  cross-sections  of 
the  same  area,  it  was  shown  that  the  anisotropic  bar  with 
maximum  torsional  rigidity  had  an  elliptic  shape, 

The  optimal  nonhomogeneity  in  torsion  for  perfectly 
plastic  bars  (yield  limit  was  the  control  function)  was 
considered  by  Mi oduchowski [ 1 3 , 1 4 ] .  He  employed  the  method 
of  incomplete  gradient  and  local  variation  on  prismatic  bars 
with  square,  circular  and  elliptic  cross-sections.  His 
study  revealed  that  'for  solutions  ensuring  maximum  load 
carrying  capacity,  zones  of  higher  yield  limit  occur  in 
regions  where  plastic  zones  would  first  occur  in  a 
homogeneous  bar  of  the  same  cross-section'  .  It  was  further 
shown  that  'the  zones  of  higher  yield  limit  form  around  the 
center  of  the  lateral  surfaces  when  the  difference  between 
two  material  constants  is  small  while  with  increasing 
difference  the  zones  change  to  form  a  belt  shape  around  the 
center  of  the  cross-section' .  These  conclusions  are 
remarkably  similar  to  those  presented  by  K losowi cz [ 9 ] . 

It  should  be  mentioned  here  that  in  the  works  of 
Klosowicz  and  M i oduchowski ,  the  control  function  (either  the 
shear  modulus  or  yield  limit  function)  is  assumed  continuous 
over  the  cross-section  such  as  might  occur  during  a 
mechanical  or  heat  treatment  of  the  bar.  It  is  however 
almost  impossible  by  this  means  to  obtain  a  desired 
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distribution  of  material  constants.  Furthermore,  it  is 
unlikely  that  the  average  shear  modulus  or  yield  limit  has 
any  significance  for  practical  engineering  problems. 

In  the  present  investigation,  the  problem  of 
optimization  of  the  same  nature  as  found  in  [2]  and  [9,10] 
is  studied.  However  the  shear  modulus  in  this  problem 
varies  in  a  jump-like  manner  which  often  occurs  in  practice 
when  a  bar  is  composed  of  prismatic  parts  of  two  completely 
different  materials. 

The  proportion  of  reinforcement  instead  of  the  average 
shear  modulus  is  prescribed  as  a  constraint  and  the 
optimizing  step  function  of  shear  modulus  which  satisfies 
all  the  state  equation,  boundary  conditions  and  the 
constraint  will  be  sought.  The  iterative  algorithm 
constructed  for  this  purpose  is  very  simple  and  easy 
compared  with  those  found  in  [2]  and  [10]. 

This  study  considers  only  the  1  inear ly-elast ic  behavior 
of  the  bars,  however,  it  is  believed  that  the  perfectly 
plastic  behavior  would  be  similar. 


2.  FORMULATION  OF  THE  PROBLEM 


2.1  Sa int -Venant  Torsion 

A  general  solid  prismatic  bar  which  is  subjected  to 
torsion  is  shown  in  Figure  2.1.  Under  the  assumptions  of 
Saint-Venant  for  bars  of  this  type, 
each  cross-section  warps 
the  same  way  and 
each  rotates  by  an 
angle  propor tonal  to 
its  di stance  from 
the  origin  with  no 
in-plane  distortion. 

This  is  expressed 
ana lyt ica 1 ly  by  Y 


z 


vx  ~  @yz'  Fig.  2.1  Coordinate  system  of 

a  prismatic  bar 

(2. 1 )  v  =  Gzx , 

vz  =  G'I'i  xty). 


where  defines  the  warping  and  is  a  function  of  x  and  y 
alone,  and  G  is  the  angle  of  twist  per  unit  length. 

The  only  nonvanishing  stress  components  are  derived  by 
the  stress-strain  law  as 


(2.2) 


Tzx=  GG^-y)' 
^zy  "  +x) . 
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A  substitution  of  these  values  into  the  equations  of 
equilibrium  yields  the  Laplace  equation  as 

g20  a20 

(2.3)  grad 20  =  ^  =  0. 

The  boundary  condition  of  the  tract  ion- free  lateral 
surface  yields 

izk-y>U-<dS+x>£s  =  o. 


while  the  shear  stresses  on  the  ends  result  in  a  twisting 
moment  as 


(2.5)  M  =  (x2+y2+x|y  -y^]dxdy. 

R 

The  Prandtl's  stress  function  0  is  defined  in  terms  of 
warping  function  or  the  shear  stresses  by 


(2.6) 


60_  _ 

ay 

_ 

ax  " 


Tzx  =  Ger^-y), 


zy 


ax 


With  this  definition, 
the  torsion  problem 
reduces  to  a 
Di r ichlet- type  problem 
of  finding  a  function 
< P(xfy )  which  satisfies 


Fig.  2.2  Shear  stresses  on  a 
cross-sect  ion 


(2.7) 


grad 2 0  =  -2GB  in  Rf 


and  the  boundary  condition  on  the  lateral  surface 
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(2.8)  0=0  on  C. 


The  expression  for  the  twisting  moment  in  terms  of  0  is  then 
(2.9)  M  =  2  0dxdy. 


It  can  be  shown  that  these  two  formulations  above  are 
derived  from  the  principles  of  minimum  potential  energy  and 
of  minimum  complementary  energy.  The  formulation  in  the 
present  study  will  be  based  on  the  latter  which  states: 
'Among  all  the  sets  of  admissible  stresses  that  satisfy  the 
equations  of  equilibrium  and  the  boundary  conditions,  the 
set  of  actual  stress  components  makes  the  total 
complementary  energy  an  absolute  minimum'  . 

The  complementary  energy  is  expressed  in  terms  of  the 
stress  function  as 


t2-10>  nc  If'2  +,%y>2] dxdy  -  2ejjj>dxdy. 

R  R 


The  solution  0  which  minimizes  7TC  is  also  the  unique 
solution  of 


ff  L  grad  0- grad  T]  dxdy  =  26  [  f  TJdxdy 
J  JrG  jjr 


(2.11 ) 


for  all  admissible  7). 

Here,  grad0  grad7/ denotes  the  scalar  product  of  two  vectors 


(2.12)  grad0-gradT)  = 


d0  dT)  .  60  BJL 
ex  ax  ay  ay 


and  the  admissible  function  7]  is  in  the  class  of  functions 
that  are  continuous  over  the  cross-section  and  satisfy  the 
boundary  conditions 


- 
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(2.13)  Tj  =  0  onC. 

It  can  be  easily  shown  that  eq . ( 2 . 1 1 )  is  equivalent  to 
the  compatibility  equation.  Consider  the  LHS  functional  in 
eq . ( 2  .  1  1  ) 


(2. 14) 

J 1  i  grad  0  •  grad  T)  dxdy 

Taking  the  eq . ( 2 . 1 3 )  into  consideration,  this  can  be 
integrated  as 


(2.  15) 

f  r  i ,  a20  a20  ^ 

J  JR  C  3* 2  ay 2  ^  dxdy . 

Substitution  of  eq . ( 2 . 7 )  into  eq . ( 2 . 1 5 )  finally  yields 


(2. 16) 

JJ  l(  206)71  dxdy  =  2ejj  7)dxdy 

R  R 

for  a  1 1 

admi  ssible  T) . 

It 

is  sometimes  more  convenient  to  use  the  eq . ( 2 . 1 1 ) 

for  the  purpose  of  theoretical  developments  and  proofs.  The 
formulation  in  the  next  section  employs  both  of  eqs . ( 2 . 1 0 ) 
and  (2.11). 


. 


2.2  Formulation  of  the  Optimization  Problem 


1 1 


Figure  2.3  illustrates 
composed  of  a  few  prismatic 
two  different  materials 
of  shear  moduli  G^  and  G 2. 
The  problem  is  to  find 
the  optimal  distribution 
of  those  components  such 
that  the  maximum 
torsional  rigidity 
can  be  achieved. 

The  procedure  fol lowed 
is  similar  to  that  of 
Klosowicz  and  Lurie! 10]. 

The  nonhomogeneity  of 
modulus  function 


a  nonhomogeneous ,  prismatic  bar 
components  but  made  of 


Fig.  2.3  Cross-section  of  a 
nonhomogeneous  bar 


bar  is  described  by  the  shear 


(2.17) 


Here,  it 


G(x,y)  =  | 
is  assumed 


G-j  in 
G  2  i  n 

that 


Region  1 
Region  2. 


(2.18)  0  <  G]  <  G2. 


The  proportion  of  rei nforcement  is  prescribed  as  a 
constraint 


(2.19) 
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where  A,  A ^  and  A 2  denote  the  areas  of  the  cross-section, 
region  1  and  region  2,  respectively. 

The  objective  is  to  achieve  the  maximum  torsional 
rigidity  or  the  maximum  twisting  moment  for  a  given  angle  of 


twist 

(2.20) 

M  =  2[\<p  dxdy. 

This  Prandtl's  stress  function  0  must  minimize  the 
complementary  energy 


( 2.21 ) 

7 Tc  ~W\  £  ( grad0) 2  dxdy  -  20 J f  <pdxdy. 

R  R 

If  the  specific  compliance  is  defined  as 


(2.22) 

u(x’y>  =G(x,y) ' 

then  eq . ( 2 . 2 1 )  can  be  rewritten  as 


(2.23) 

/(c  =2//  u  2dxdy  -  20 jf  (pdxdy. 

The  stress  function  in  eq.(2.23)  is  the  unique  solution 
of 


(2.24) 

JJ  u  grad 0-  grad  Tf  dxdy  =  2oJJr)dxdy 

R  R 

for  a  1 1  admi  ssible  T] 
If  0  =  V  then 


(2.25) 

iff  u  (grad<p)2dxdy  -  20 ff  (Pdxdy 

<6/ Jr  'Jr 

=  -off  0 dxdy 

JJR 

Hence  the  problem  can  be  stated  as 


(2.26) 


13 


max.  off  (pdxdy. 


uf(p 


where  u  and  satisfy  eq.(2.24)  with  the  constraints 


(2.27)  u(x,y)  = 


U]  in  Region  1 
U2  in  Region  2 


where 


(2.28)  0  <  U2  ^  U] 


3.  SOLUTION  PROCEDURE 


3.1  Necessary  Conditions  for  Optimality 

The  solution  for  the  problem  stated  in  the  previous 
chapter  does  not  appear  possible  by  analytical  means. 
However,  since  the  objective  of  the  present  study  is  to 
develop  a  simple  and  versatile  method  which  is  applicable  to 
any  shape  of  cross-sections,  the  F.E.M.  was  used  in 
conjunction  with  an  iterative  procedure. 

Although  no  formal  proof  of  sufficient  and  necessary 
conditions  as  well  as  the  existence  and  convergence  of  the 
solution  has  been  developed,  these  are  accomplished  in  the 
practical  manner  in  the  next  chapters. 

In  order  to  develop  an  iterative  algorithm  consider  the 
functions  ur  and  0r  are  computed  and  ur+1  and  0r+1  in  the 
next  iteration  are  to  be  computed. 

If 


(3.  1 ) 


ur»i  =  ur  +  Auc 
0r*l  =  0r  +  A(t> r  - 


then  the  objective  functional  for  the  r+1-th  step  becomes 


(3.2)  Jr+1  =  JJ^(ur+Aur){grad((pr  +  A<pr)}2dxdy 

=  26  J J^(Qr+A(pr)dxdy . 


But  since 
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. 
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(3.3>  J f  (ur+Aur  )grad(<Pr+A0r).  gradT)  dxdy 
^  *  R 

=  20  f  [  Tjdxdy 
JJR 


=  J J  ux  g rad grad  T)  dxdy , 

R 

for  all  admissible  V, 

the  functional  in  eq . ( 3 . 2 )  can  be  expressed  as 

(3.4)  Jr  +  1  =  JJ  urgrad0r- grad(0r +A0r  )dxdy 

R 

=  3r+JJ  urgrad0r-gradA0rdxdy . 

R 

Again  from  the  identity  eq.(2.24) 

(3.51  J J  L/r+1  grad(0r+A0r)  ■  grad0rdxdy 

R 

=  20  Jf  0rdxdy 

R 

=  ffRUr  (grad0r  )2dxdy , 
the  relationship  is  obtained  as 
(3.6)  JJ  ur  grad 0rgradA0r  dxdy 

+  J J^Aurgrad(0r+A0r )  •  grad0rdxdy 


=  0. 


By  substitution  of  eq.(3.6)  into  eq . (3.4) ,  the 
objective  functional  finally  becomes 

=  Jr  -  JJ  Aurgrad(<t>r  + A0r) -grad0rdxdy 
R 


(3.7) 


Therefore,  the  increase  of  the  value  of  the  functional  can 
be  expressed  as 
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(3.8)  Adr  =  -  jj  Aur( grad <Pr)2dxdy 

-  jj  Aurgrad(pr  gradA(prdxdy. 


If  grad^0ris  negligible  compared  with  grad0r,  Adr  can 
be  approximated  as 


(3.9) 


ur(grad  <J>r)  2dxdy. 


The  necessary  conditions  for  achieving  the  maximum 
torsional  rigidity  is  that  AJr  in  eq . ( 3 . 9 )  should  be 
positive,  or 

(3.10)  jj  Au f( grad 0r ) 2dxdy  <  0. 

R 

Consider  a  cross-section  of  a  bar  divided  into  a  number 
of  constant  strain  triangles  -  see  Figure  3.1.  And  let 
t/r^1  be  different  from  ur  in  the  elements  /  and  j  alone, 
i  .  e .  , 


(3.11 ) 


u2~u  1 
uru  2 


0 


in  element  / 
in  element  j 
e 1 sewhere . 


The  integration  in  eq . ( 3 . 1 0 )  can  now  be  performed  as 
(3.12)  Adr  s  ( U]  -u2 )  [(grad0r)2.  ~(gradQr)2.  } 

Since  u 2  <  (or  <  G2)  ,  Adr  is  positive  when 


(3.13)  (grad(Pr)2.  >  ( grad(pr) 2  ■ 


It  can  be  noted  here 


dd)  dd> 

(3.14)  ( grad<p) 2  =  > 2  +  (p ) 2 


=  T  2  +T  2 
'zy  'zx 


Fig.  3.1  Interchange  of 

elements  /  and  j 


Hence  the  conclusion 
may  be  stated  as: 

In  order  to  improve  the  torsional  rigidity  of  a 
certain  nonhomogeneous  bar,  it  is  necessary  and 
sufficient  in  the  approximate  sense  to  change  the 
shear  modulus  distribution  such  that  the  higher 
stress  area  in  the  more  flexible  material  region 
should  be  interchanged  with  the  lower  stress  area 
in  the  stiffer  material  region. 
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3.2  Iterative  Algorithm 

To  determine  the  maximum  torsional  rigidity,  an 
iterative  algorithm  is  constructed  based  on  the  conditions 
described  previously  as  follows: 

i)  Select  an  initial  shear  modulus  distribution  such 
that  it  satisfies  the  constraint,  eq.(2.19), 

II)  Compute  the  stresses  in  each  element  as  well  as 
the  torsional  rigidity  of  the  cross-section, 

iii)  Select  two  elements;  one  of  the  highest  stress 
from  region  1,  the  other  of  the  lowest  stress 
from  region  2, 

iv)  Define  the  new  shear  modulus  function  Gr*i  which 
is  locally  different  from  Gr  by  interchanging  the 
shear  modulus  in  .those  two  elements, 

v)  Repeat  the  steps  ii),  iii)  and  iv)  until  the 
maximum  torsional  rigidity  is  obtained. 

The  initial  shear  modulus  function  can  be 
arbitrary  as  long  as  it  satisfies  the  constraint. 
From  then  on  the  proportion  of  rei nforcement  will 
be  Kept  the  same  by  interchanging  a  pair  of 
elements  of  the  same  area. 

When  the  hybrid  stress  formulation  is  used 
in  step  ii)  as  in  the  case  of  the  present  study, 
the  shear  stresses  in  each  element  are  calculated 
from  the  values  of  warping  function  instead  of 
the  Prandtl's  stress  function. 


In  step  iii),  the  number  of  elements  to  be 
interchanged  needs  not  be  two,  as  several  pairs 
of  elements  can  be  selected.  In  the  present 
study,  five  pairs  of  elements  were  selected  in 
the  first  step.  This  number  was  then  reduced  as 
the  solution  approached  the  optimum. 

The  interchange  of  the  shear  moduli  of  the 
elements  in  step  iv)  causes  modification  of  the 
equations  cor respondi ng  to  the  nodal  points 
around  these  elements.  This  modified  set  of 
equations  is  then  solved  in  step  ii) 

The  iteration  stops  when  either  no  element 
is  chosen  in  step  iii)  or  the  interchange  of  even 
two  elements  results  in  a  decrease  of  the 
torsional  rigidity. 


4.  FINITE  ELEMENT  AND  COMPUTER  PROGRAM 


4.1  Hybrid  Stress  Formulation 

The  hybrid  stress  approach  suggested  by  Yamada  et 
al.[22]  is  employed  for  the  finite  element  torsion  analysis 
in  this  study. 

The  cross-section  of  a  prismatic  bar  is  divided  into  a 
number  of  constant  strain  triangular  elements,  and  these 
elements  are  classified  into  the  'inner'  and  'boundary' 
elements.  A  boundary  element  is  an  element  which  possesses 
a  side  along  the  boundary  of  the  cross-section.  Figure  4.1 
shows  a  typical  element  division  and  a  boundary  element. 

The  nodal  points  are  numbered  counter-clockwise  i-j-k 
so  that  a  boundary  element  should  have  j  and  k  nodes  on  the 
boundary. 

The  stress  function  is  assumed  inside  an  element  and  is 
used  for  the  stress  boundary  condition,  while  the  warping 
function  is  defined  along  the  three  sides  of  an  element  so 
that  the  displacement  is  naturally  continuous  over  the 
cross-section.  Both  the  stress  and  warping  are  assumed 
linear  in  the  Constant  Strain  Triangular  elements  (Figure 
4.2)  . 

In  this  case  of  hybrid  stress  approach,  the  relation 
between  the  stress  and  the  warping  function  is  then  derived 
from  the  principle  of  minimum  complementary  energy.  The 
complementary  energy  expression  for  a  triangular  element 
prism  with  unit  thickness  is 


20 
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(4. 1 )  llc=  U  -  W. 


BOUNDARY  ELEMENT 


© 


(a)  Element  classification  (b)  A  boundary  element 
Fig.  4.1  'Inner7  and  'Boundary'  elements 


The  first  term  in  the  expression  of  the  complementary 
energy,  eq.(4/1)  is  the  strain  energy  stored  in  an  element 
prism  with  unit  thickness 

(4.2)  U  =  IjJJ  (0}T  [D]{Q}dxdy 

A 

where  {0}  is  defined  as 


(4.3a) 


for  an  inner  element.  Here  in  eq.(4.3a)  A  denotes  the  area 
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of  the  element  and  a|=X|<-Xj,  b-{  =y-  -yk  in  cyclic  order  of  i  ,j 
and  K.  For  a  boundary  element  with  a  side  j-k  along  the 
boundary,  0-  =  0k=  0  and  eq.(4.3a)  simplifies  to 


or  in  P i an' s [ 1 7 ]  original  notation 


(4.3c)  {Q}  =  [P]{{]}. 

In  eq . ( 4 . 3c ) , 


(4.4a) 


(4.4b) 


IP) 

[P] 


an  inner  element, 


an  for  a  boundary  element 


and 


(4.5a) 

<0}  =  < 

K-nl 

[^r0kj 

for  an  inner  element, 

(4.5b) 

<fi}  =  <t>: 

for  a  boundary  element. 

The  matrix  [D]  in  eq . ( 4 . 2 )  denotes  the  inverse  of  the 
shear  modulus 


(4.6) 


ID]  =  k 


1  0 
0  1 


y 


s 


Fig. 


4.2  Assumed  warping 


f unct i on 
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If  the  integration  in  eq.(4.2)  is  carried  out,  the 
strain  energy  term  becomes 

(4.7)  U  =  ^jj  {CJ}T  [D]iO}dxdy 

A 

where  [H]  is  a  matrix  defined  by 

(4.8)  [H]  =  /I [ P]T  ID]  [P]  . 

The  second  part  of  the  complementary  energy  is  the  work 
of  the  surface  tractions  over  the  boundary  surfaces  of  an 
element  prism  where  the  displacement  is  defined 

(4.9)  VI  =  j{T)U^b}ds 

where  { T } T  =  <M  Tkj  Tjj  >  is  the  vector  of  twisting 

moment  and  normal  components  of  shear  stresses  on  three 
sides  of  an  element;  ^ki  ^ij  >  the  vector 

of  corresponding  angle  of  twist  and  warping  asssociated  with 
the  sides  of  the  element.  The  work  of  the  external  loads  W 
can  now  be  integrated  to 

(4.10)  iV  =  {(1}T  [P]T  [  i. ]  . 

Here , 

(4.11)  {q^} T  =  <0  '/'j  ^k> 

and 
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b\  bj  bk 
a\  a j  ak  .  ' 

The  complementary  energy  is  then  given  from  eqs . ( 4 . 1 ) , 
(4.7)  and  (4.10)  as 

(4.13)  7TC  =  ^  {f3>T  tW]  {(3>  -  {B}T[P]TUHq^} 


(4.12) 


ui  -  i 


-2y0A 


2x0A 


Differentiating  7TC  with  respect  to  {HI,  and  equating 
the  result  to  zero,  the  realation  is  established  as 

(4.14)  {ft}  =  [H]"1  [P]T  [  /-  ]  {q^}  . 

The  strain  energy  U  can  be  expressed,  by  definition,  in 
terms  of  stiffness  matrix  as 


(4.15)  U= 

where  [k]  is  the  stiffness  matrix. 

From  the  eqs. (4.7)  and  (4.14),  and  by  noting  that  the 
martix  [H]  is  symmetric,  the  alternative  expression  for  U 
can  be  derived  as 

(4.16)  U  =  |{B)T  [HUB} 

=  l  <VT  l*-]1  [P]  [wr'  I P]JU]{q^}. 

The  element  stiffness  matrix  associated  with  the  hybrid 
approach  is  now  obtained  by  equating  the  above  two 
expressions,  eqs. (4. 15)  and  (4.16)  for  U  as 


■* 


26 


(4.17a)  [k]  =  [L]J  (P}[H}-]  [P)J(L] 


G_ 

4A 


4A2(x2+y2  ) 
o  Jo 

2A(a.x  -b.y  )  a2+b2 

10  10  II 

2A(  a-x-b-y  )  a-a-+b-b.  a2+b 2 
jo  jo  ‘  J  1  J  J  J 


Sym. 


24iakxo-Vo ;  aiak+bibk  a,ak+bjbk  ak+bd 


for  the  inner  elements,  and 


(4. 17b) 


[A]  =  TT 


GA 

Jk 


(a  y  +b-  x  ) 2 
1*0  10 

0  0 

a-yn+b-xn  0 

jo  J  0 

~akyo'bkxo  0 


1 

-1 


Sym. 


for  the  boundary  elements.  Here,  /•£  denotes  the  squared 
length  of  the  j-k  side. 

Furthermore,  the  stiffness  matrix  [k] ,  by  its 
definition,  connects  the  generalized  force  and  displacement 
vectors  as 


(4.18a) 


> 

M 

/  \ 
0 

F: 

1 

=  1*1  <  > 

F, 

/k, 

{  / 

or 


(4. 18b) 


{<?}  =  IfcHq*} 


27 


where  F j  ,  F j  and  represent  the  resultant  shear  forces  on 
the  sides  j-k,  k- i  and  i-j,  respectively.  When  the 
generalized  force  vectors  of  all  the  elements  are  assembled 
into  the  global  vector,  the  scalar  sum  of  these  forces 
around  a  nodal  point  should  be  equated  to  zero  so  that  the 
requirement  of  equilibrium  may  be  fulfilled.  In  other 
words,  eq.(4.18a)  for  the  whole  cross-section  is  assembled 
as 


(4. 19) 


The  vector  {Fn}  is  a  zero  vector 


(4.20)  {Fn }  =  {0} ,  n =1,2,3 


f  •  •  f 


number  of  nodes. 


This  means 


(4.21)  ikmf}6  +  [k ff  ]{^n}  =  {0} 


so  that  the  torsional  rigidity  can  be  obtained  as 


k 


mm 


{kmf}T  [k  ff  ]  1  {^mf} 


The  stresses  in  each  element  can  be  calculated  from 


9 


(4.23) 
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Fig. 


Fig.  4 


1.3  Shear  forces  on  the  three  sides  of  an  element 


(D  © 


© 


[H  ELEMENT 


4  Local  and  global  numbering  system  of  nodes  and 
elements 
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4.2  Flow  Charts  of  the  Computer  Program 

The  computer  program  for  the  present  study  of 
optimization  is  written  in  FORTRAN,  and  consists  of  MAIN  and 
four  major  subroutines,  namely;  DATA,  ASSEMB ,  SOLVE  and 
STRESS.  Flow  charts  for  each  of  these  subroutines  are 
included  in  the  next  few  pages.  The  complete  program  is 
attached  in  the  appendix. 
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i)  MAIN  controls  the  sequence  of  computation  by  calling  the 
subroutines,  and  makes  the  decision  whether  or  not  to  carry 
out  further  iteration  and  if  necessary,  how  many  elements 
are  to  be  interchanged. 


' 
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ii)  Subroutine  DATA  reads  the  given  nodal  and  element  data, 
generates  and  completes  the  data  for  intermediate  nodes  and 
elements,  prints  out  complete  nodal  and  element  data,  and 
determines  the  number  and  band  width  of  the  equations  to  be 
solved . 
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iii)  Subroutine  ASSEMB  takes  an  element  in  turn,  constructs 
an  element  stiffness  matrix,  assembles  the  global  matrix  and 
then  applies  the  displacement  boundary  condition. 
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iv)  Subroutine  SOLVE  solves  eq.(4.23)  for  the  values  of 
warping  function  at  the  nodal  points  by  using  the  Gaussian 
elimination  technique  then  finds  the  torsional  rigidity  for 
the  present  shear  modulus  distribution. 
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v)  Subroutine  STRESS  calculates  the  shear  stresses  in  each 
element,  selects  the  elements  to  be  interchanged,  and 
modifies  the  stiffness  matrix  according  to  the  new  shear 
modulus  distribution. 


(  start) 


Compute  the  nondimens iona 1  shear 
stresses  TXZ,  TYZ  and  their  square 
sum  SQT  in  each  element 


Select  an  element  /  with  the 
largest  SQT  in  weak  material  zone 


Select  an  element  j  with  the 
samllest  SQT  in  stiffer  material  zone 


(return) 

- ft - 


No 


Interchange  materials 
in  element  /  and  j 


Modify  the 
stiffness  matrix 


if  No. 
of  interchange 
<NECH 


Yes 


5.  NUMERICAL  TESTS  AND  RESULTS 


5.1  Testing  of  the  Numerical  Procedure 

As  the  solution  procedure  outlined  in  chapter  3  is  not 
completely  rigorous,  it  is  necessary  to  demonstrate  the 
uniqueness  and  convergence  of  the  solution  as  well  as  the 
reliability  of  the  computer  program  by  considering  special 
test  cases.  In  the  following  subsections  the  results  of  the 
F.E.M.  analysis  for  a  bar  with  square  cross-section  are 
presented  and  compared  with  Known  solutions. 

5.1.1  The  case  of  a  homogeneous  bar 

As  a  first  numerical  example,  the  finite  element 
program  based  on  the  hybrid  stress  approach  was  applied  to  a 
homogeneous  square  bar.  Due  to  the  symmetry  of  the  square 
cross-section,  only  one-eighth  of  it  was  divided  into  one 
hundred  elements  and  sixty-six  nodes.  Figure  5.1  shows  the 
element  division  of  a  square  bar. 

The  calculated  torsional  rigidity  for  a  square  bar  of 
side  7  in  non-dimensional  form  is 
J  =  g^j7  =  0.1404 

while  the  analytical  solution  is  given  in  [20]  as 
J  =  =  0.1406 

The  torsional  rigidity  obtained  with  different  numbers 
of  element  divisions  is  shown  in  Table  5.1  and  Figure  5.3. 
Figure  5.3  also  shows  the  torsional  rigidity  calculated  by 
the  displacement  approach. 
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Fig.  5.1  Element  division 
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Fig.  5.2  Contours  of  constant  warping  on  the  cross-section 
of  a  homogeneous  square  bar 
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Fig.  5.3  Torsional  rigidity  of  a  homogeneous  square  bar 
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In  addition,  the  contours  of  constant  warping  function 
are  plotted  in  Figure  5.2. 


Table  5.1  Torsional  rigidity  of  a  homogeneous  bar 


No.  of  -j2  2 2  3^  42 

elements 


5 2 


Ana lyt ica 1 
solution 


J  0.111  0.128  0.135  0.138  0.139  0.140  0.140  0.141 


5.1.2  The  case  of  a  simple  nonhomogeneous  bar 

The  results  of  the  computer  program  in  the  case  of 
nonhomogeneous  bars  were  also  checked  for  a  square  bar 
composed  of  two  equal  parts  of  different  materials.  In  this 
case,  one  half  of  the  square  had  to  be  taken  into 
consideration  -  see  Figure  5.4.  The  calculations  were  done 
with  one  hundred  elements  and  sixty-six  nodes  and  shear 
modulus  ratios  of  1.0,  1.5,  2.0,  3.0,  5.0  and  10.0.  The 
results  are  tabulated  and  plotted  in  Table  5.2  and  Figure 
5.4.  The  analytical  solutions  calculated  from  a  formula  in 
[16]  are  also  presented  for  the  purpose  of  comparison. 

Table  5.2  Torsional  rigidity  of  a  simple  nonhomogeneous  bar 


G 1 

1  .0 

1.5 

2.0 

3.0 

5.0 

10.0 

FEM  Result 

0.139 

0.170 

0.  195 

0.236 

0.306 

0.458 

0.141 

0.172 

0.  197 

0.239 

0.311 

0.466 

G,Ol* 
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Fig.  5.4  Element  division  and  the  torsional  rigidity  of  a 
simple  nonhomogeneous  square  bar 
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5.1.3  The  convergence  and  uniqueness 

To  test  for  the  convergence  to  the  optimal  solution  as 
well  as  its  uniqueness,  a  square  bar  composed  of  75%  G  and 
25%  G  materials  with  shear  modulus  ratio  of  1.5  was 
considered.  An  eighth  of  the  cross-section  was  again 
divided  into  one  hundred  elements  and  sixty-six  nodes. 
Starting  from  completely  different  initial  distributions, 
the  optimal  distribution  was  computed.  Figure  5.5  and  5.6 
show  the  sequence  of  shear  modulus  distributions  which  were 
calculated  in  each  iteration,  and  the  torsional  rigidity  for 
each  iteration  is  shown  in  Figure  5.7.  Note  that  five 
elements  were  interchanged  in  the  first  step,  the  computer 
program  then  reduced  this  number  by  one  whenever  the 
increment  of  torsional  rigidity  became  negative. 
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Fig.  5.5  Sequence  of  iteration  starting  with  stiffer 
material  in  elements  No. 1-25 
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Fig.  5.6  Sequence  of  iteration  starting  with  stiffer 
material  in  elements  No. 76-100 


vie'D 
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Fig.  5.7  Calculated  torsional  rigidity  for  different  initial 
shear  modulus  distributions 
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5.2  Solutions  for  optimal  nonhomogeneity 

Optimal  solutions  for  the  square,  nonhomogeneous  bars 
are  considered  for  various  shear  modulus  ratios,  and  the 
comparisons  of  these  optimal  solutions  to  those  of  a 
circular  bar  are  made, 
i )  Solutions  for  a  square  bar 

A  square  bar  of  side  7  and  unit  length  was  considered 
for  various  shear  modulus  ratios  and  proportions  of 
rei nforcement .  Figure  5.9(a),  (b)  and  (c)  show  the 
optimizing  distribution  for  three  different  proportions  of 
reinforcement,  0.25,  0.5  and  0.75  each  for  the  shear  modulus 
ratios  of  1.5,  2.0,  3.0,  5.0,  10.0  and  20.0  The  maximum 
torsional  rigidity  obtained  in  each  case  is  tabulated  in 
Table  5.3  and  plotted  in  Figure  5.8. 

Table  5.3  Maximum  torsional  rigidity 


G  i 

1.5 

2.0 

o 

CO 

5.0 

10.0 

20.0 

0.25 

0.174 

0.205 

0.265 

0.386 

0.683 

1  .238 

p 

0.5 

0.  193 

0.245 

0.348 

0.556 

1  .072 

2.  102 

0.75 

0.206 

0.271 

0.402 

0.663 

1  .316 

2.621 

. 


"  Gi0l4 
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Fig.  5.8  Optimal  nonhomogeneity  of  a  square  bar 
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Fig.  5.9  The  maximum  torsional  rigidity  for  various  shear 
modulus  ratios  and  proportions  of  reinforcement 


(a)  p  =  0.25 
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(b)  p  =  0.5 
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(c)  p  =  0.75 
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5.3  Optimal  Solutions  of  a  Circular  Bar 

To  compare  the  optimal  solution  of  the  square  bar  with 
that  of  a  circular  bar,  a  circular  bar  with  radius  a  which 
has  the  same  cross-sectional  area  is  considered. 

Since 

(5.1 )  TCa2  -- 
then 


(5.2) 


a 


Fig.  5.10  Optimum  nonhomogeneous  circular  bar 


It  can  be  easily  shown  that  the  optimal  nonhomogeneous 
circular  bar  is  of  the  form  as  illustrated  in  Figure  5.10, 
of  which  the  torsional  rigidity  is  given  as 


S-  ( 7^(aA-b* )  +  bA ) . 


The  proportion  of  rei nforcement  is  simply  the  ratio  of 
the  ring  area  to  that  of  entire  cross-section 


(5.4) 


r. 
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The  substitution  of  eq.(5.4)  into  eq . (5.3)  gives  the 
expression  of  the  maximum  torsional  rigidity  in  terms  of  p 
as 

(5.5)  J  =  +  (1-p)2}. 


For  a  fixed  shear  modulus  ratio,  the  optimal  torsional 
rigidity  is  a  quadratic  function  of  p  and  is  a  linear 
function  of  shear  modulus  ratio  for  a  given  p. 

Now  define  the  efficiency  of  optimization  as 


(5.6) 


eff . 


<^op 

p(d2~dy  ) 


where  d  ,  d  and  d  denote  the  torsional  rigidities  of 
homogeneous  bars  made  of  G^ ,  G ^  materials  and  of  optimal 
nonhomogeneous  bar,  respectively.  This  efficiency  of 
optimization  can  be  understood  as  the  ratio  of  the 
improvement  in  the  torsional  rigidity  to  the  expected 
improvement  with  the  same  amount  of  stiffer  material  if  the 
increment  of  the  torsional  rigidity  were  proportional  to  the 
amounts  of  stiffer  material  being  reinforced.  From  the 
eqs . ( 5 . 3 )  ,  (5.5)  and  (5.6),  the  expression  of  the  efficiency 
of  optimization  can  be  derived  in  terms  of  p  for  a  circular 
bar  as 


(5.7)  eff.  =  2-p. 

These  efficiencies  for  a  circular  and  a  square  bar 
depending  on  the  proportion  of  rei nforcement  can  be  seen  in 
F igure  5.11. 


EFFICIENCY,  eff. 


52 


PROPORTION  OF  REINFORCEMENT,  P 


Fig.  5.11  Efficiency  of  optimization 
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1.  DISCUSSION  AND  CONCLUDING  REMARKS 


1.1  Discussion  of  the  Solution  Procedure 

The  F.E.M.  solution  procedure  used  in  this  study  is  a 
compromise  between  the  using  a  very  simple  element  combined 
with  a  trial  and  error  procedure  for  optimization,  and  that 
of  using  a  complex  element  combined  with  a  very  fine  element 
mesh.  The  simplest  approach  may  be  too  crude  to  find  an 
answer  close  to  the  actual  optimum,  while  the  latter 
approach  would  be  quite  expensive  to  use. 

As  seen  in  Figure  5.3,  the  calculated  torsional 
rigidity  for  a  homogeneous  square  bar  converges  to  the 
analytical  solution  with  increasing  element  numbers.  The 
torsional  rigidity  obtained  by  the  displacement  approach 
instead  of  the  hybrid  stress  approach  is  also  shown  in 
Figure  5.3.  It  is  interesting  to  note  that  the  torsonal 
rigidity  by  the  displacement  approach  converges  from  above 
the  Known  solution  (upper  bound),  while  the  torsional 
rigidity  by  the  hybrid  stress  approach  converges  more 
rapidly  from  below  (lower  bound).  Furthermore,  the  stress 
solutions  by  the  hybrid  stress  approach  should  also  be  more 
accurate  than  those  by  the  displacement  approach.  Since  the 
stress  in  each  element  is  to  be  compared  in  order  to  select 
a  pair  of  elements  for  interchange,  the  hybrid  stress 
approach  was  employed  in  the  present  study.  A  more  detailed 
discussion  of  the  use  of  displacement  and  hybrid  stress 
approach  can  be  found  in  [3,22]. 
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In  the  actual  computer  program,  the  shear  stresses  are 
nond i mens iona 1 i zed  and  then  compared,  i .e. ,  the  the  shear 
stress  divided  by  the  shear  modulus  is  used  for  selecting 
elements  to  be  interchanged.  This  is  mainly  for  the  reason 
that  the  elements  in  the  stiffer  material  region  can  always 
have  higher  stresses  when  the  shear  modulus  ratio  is  high 
enough  so  that  the  interchange  of  the  shear  modulus 
distribution  can  not  be  accomplished.  Also  by  interchanging 
more  elements,  more  shear  modulus  distributions  may  be 
considered  with  the  result  that  the  chance  of  missing  the 
best  one  may  be  reduced. 

Five  pairs  of  elements  at  the  most  are  to  be  chosen  in 
the  first  iteration  and  this  number  is  reduced  by  one 
whenever  the  increment  of  torsional  rigidity  becomes 
negative.  This  number  can  be  arbitrarily  set,  in  fact,  one 
may  choose  any  number  of  element  pairs  to  begin.  However, 
if  too  many  elements  are  interchanged,  the  assumption  that 
the  change  of  stresss  function  is  small  may  be  violated;  if 
too  few,  it  may  cause  extremely  small  changes  in  the 
torsional  rigidity. 

The  size  of  element  in  the  present  study  was  taken  as 
constant  for  all  elements.  However,  the  computer  program 
can  be  modified  so  that  it  can  take  different  size  elements. 
It  would  be  possible  therefore  to  improve  the  accuracy  of 
solutions  by  having  a  finer  element  mesh  in  the  higher 
stress  region  around  the  center  of  the  lateral  surface  of 


the  bar . 
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6.2  Discussion  of  Results 

The  test  results  of  a  homogeneous  square  bar  show  the 
error  of  the  calculated  torsional  rigidity  by  the  hybrid 
stress  approach  being  1.2%  and  0.2%  with  25  and  100  element 
divisions,  respectively  (Table  5.1  and  Figure  5.2),  while 
the  errors  for  the  nonhomogeneous  bar  vary  from  1.2%  to  1.7% 
depending  on  the  shear  modulus  ratio  (Table  5.2  and  Figure 
5.4).  Although  half  of  the  square  instead  of  an  eighth 
section  is  taken  and  divided  into  100  elements  in  the  case 
of  a  simple  nonhomogeneous  bar,  the  size  of  the  element  is 
same  as  for  an  eighth  of  a  homogeneous  square  with  25 
elements.  It  is  therefore  expected  that  the  errors  should 
be  similar  in  these  two  cases  in  the  range  of  low  shear 
modulus  ratios.  The  error  increases  as  the  difference  of 
material  constants  increases.  This  means  that  the  results 
in  this  study  are  more  reliable  for  lower  shear  modulus 
ratios. 

Figure  5.5  and  5.6  show  that  the  optimal  solution  is 
essentially  unique  regardless  of  initial  shear  modulus 
distribution,  and  that  the  sequence  of  G  or  u  converges 
to  this  solution.  The  final  solution  in  Figure  5.5  and  5.6 
appear  to  be  different,  however  the  difference  in  torsional 
rigidity  is  very  small  (less  than  0.005%).  The  true  optimal 
solution  must  be  near  these  two  shapes  and  any  solutions 
close  to  these  two  may  be  considered  optimal  for  practical 


purposes . 
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The  approximate  uniqueness  of  solution  can  be  shown  for 
the  shear  modulus  ratios  other  than  1.5,  as  well.  However, 
as  the  shear  modulus  ratio  increases,  not  only  does  the 
accuracy  of  the  solution  decrease  but  also  it  requires  more 
iteration  steps.  Eventually,  it  becomes  almost  impossible 
to  arrive  near  the  optimal  solution  when  the  given  initial 
shear  modulus  distribution  is  too  far  from  the  final 
solution.  For  this  reason  the  example  given  is  presented 
for  a  shear  modulus  ratio  of  1.5. 

Figure  5.9  illustrates  the  obtained  maximum  torsional 
rigidity  depending  on  the  proportion  of  rei nforcement  as 
well  as  on  the  shear  modulus  ratio.  The  maximum  torsional 
rigidity  in  the  case  of  a  circular  bar  is  linearly 
proportional  to  the  shear  modulus  ratio  as  shown  in  the 
eq . ( 5 . 7  )  .  In  the  case  of  a  square  bar,  the  trend  is  almost 
the  same,i.e.,  linearly  dependent  on  the  shear  modulus 
ratio,  except  for  p=0.25.  It  is  also  seen  from  Figure  5.11 
that  the  maximum  torsional  rigidity  for  the  high  shear 
modulus  ratio  and  low  proportion  of  rei nforcement  is  less 
than  anticipated  thus  the  efficiency  of  optimization  of  a 
square  bar  falls  below  that  of  a  circular  bar.  This  is 
mostly  due  to  the  crude  shape  of  the  material  boundaries 
when  approximated  by  the  element  triangles.  As  seen  in 
Figure  5.8(a)  the  shape  of  the  rigid  rei nforcement  of  the 
optimal  solution  is  quite  smooth  in  the  lower  shear  modulus 
ratio  (1.5,  2.0)  range,  but  it  becomes  very  rough  in  the 
higher  shear  modulus  ratio  range. 
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In  order  to  investigate  this  effect,  the  optimal  shear 
modulus  distribution  for  p=0.25  and  shear  modulus  ratio  10.0 
was  considered.  Since  the  F.E.M.  solution  is  only  an 
approximate  representation  of  the  optimal  shear  modulus 
ratio,  the  actual  shape  of  the  optimal  rigid  rei nforcement 
should  be  bounded  by  a  smoooth  curve  close  to  the  saw-tooth 
shape.  If  the  solution  represented  by  the  solid  lines  in 
Figure  6.1  is  considered,  the  torsional  rigidity  is  found  to 
be 

J  =  0.696 

which  is  almost  2%  higher  than  the  previously  obtained  value 
of  0.683.  The  efficiency  of  optimization  can  then  be 
calculated  to  be 

eff.  =  1.754 

which  would  bring  the  curve  of  shear  modulus  ratio  10.0  in 
Figure  5.11  just  above  the  dotted  straight  line  for  a 
circular  bar. 

As  is  clear  from  the  above  discussion,  the  results  of 
optimal  solutions  presented  in  Chapter  5  are  approximate  and 
thus  can  be  used  only  for  a  qualitative  study  of  the 
behaviour  of  such  optimal  solutions. 

The  optimal  shear  modulus  distribution  obtained  for  the 
different  shear  modulus  ratios  and  proportion  of 
rei nforcement  are  illustrated  in  Figure  5.9(a),  (b)  and  (c). 
One  may  notice  that  the  zones  of  rigid  reinforcement  for 
lower  shear  modulus  ratios  and  lower  proportions  of 
rei nforcement  start  to  build  up  around  the  center  of  lateral 
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surface  then  with  increasing  amounts  of  stiffer  material, 
the  zones  expand  to  form  a  belt  shape  around  the  center  of 
the  cross-section.  This  observation  is  quite  consistent 
with  the  conclusions  found  in  [10]  and  [14],  although  the 
shear  modulus  or  the  yield  limit  in  their  study  is  a 
continuous  function,  not  like  the  shear  modulus  function  in 
this  study. 

The  efficiency  of  optimization  for  a  circular  bar  is 
represented  as  a  straight  line  in  Figure  5.11.  It  is  very 
interesting  that  the  efficiency  of  optimization  for  a  square 
bar  is  much  higher  than  that  of  a  circular  bar  especially  in 
the  range  of  low  shear  modulus  ratio  and  low  proportion  of 
reinforcement.  It  is  no  wonder  that  the  former  approaches 
to  latter  with  the  increasing  shear  modulus  ratio  and 
proportion  of  rei nforcement ,  because  the  zones  of  stiffer 
material  in  the  optimal  solutions  for  a  square  bar  resemble 
the  ring  shape  like  the  case  of  a  circular  bar  in  the  high 
shear  modulus  ratio  and  proportion  of  rei nforcement  range  as 
previously  discussed. 
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Fig.  6.1  Boundaries  of  the  zones  of  optimal  rigid 
reinforcement  are  made  smooth 
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6.3  Concluding  Remarks 

The  study  of  optimal  nonhomogeneity  by  K losowi cz [ 9 , 1 0 ] 
for  elastic  torsion  and  by  Mioduchowski  [13,14]  for  plastic 
torsion  revealed  interesting  characteristics  of  such  optimal 
solutions.  Zones  of  rigid  rei nforcement  start  to  appear 
around  the  center  of  lateral  surfaces  where  the  shear  stress 
would  be  the  highest  if  the  bar  were  homogeneous,  then  the 
configuration  changes  to  form  a  ring  around  the  center  of 
the  cross-section  as  the  ratio  of  material  constants  or  the 
amount  of  the  stiffer  material  increases.  It  is  more  likely 
that  the  practical  problem  often  encountered  is  of  the  type 
where  the  material  constants  vary  in  a  jump  like  manner 
rather  than  continuously.  Even  though  the  shear  modulus 
function  in  this  study  is  restricted  to  assume  either  of  two 
shear  moduli  G]  or  G 2  ,  the  optimal  distribution  for  a 
square  bar  behaves  similarly  and  hence  the  same  conclusion 
is  drawn. 

Another  aspect  of  the  optimal  solution  was  also  noted 
in  the  comparison  of  square  and  round  cross-sections.  It 
was  shown  that  a  square  bar  subjected  to  torsion  can  be  more 
efficiently  optimized;  than  a  circular  bar;  more  improvenent 
in  torsional  rigidity  with  a  smaller  amount  of  stiffer 
mater i a  1  . 

Although  only  bars  with  square  cross-section  were 
considered  in  the  present  study,  it  could  be  extended  to  the 
bars  with  cross-section  such  as  an  ellipse  or  triangles,  by 
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using  essentially  the  same  computer  program.  In  addition, 
different  sized  elements  could  be  easily  included  in  the 
procedure . 

It  is  further  believed  that  the  method  used  and  the 
solution  procedure  followed  in  the  present  study  can  be 
applied  to  the  similar  optimization  problems  in  other  fields 
such  as  ground  water  seepage  or  heat  transfer. 
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APPENDIX  -  COMPUTER  PROGRAM 

Q  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^fc  5^C  jjc  3)C  5^C  5$C  ^  ^  ^  ^  'I'  *4*  ^  ^fc  if*  'I'  ^  ^  ^  ^  41*  ^  ^  ^  *4^  4*  ^  ^  4* 

c 

COMMON  TFF( 150,30) ,EFF(3,3, 150) ,EMF(3, 150) ,NP(3, 150) , 

1  A ( 3 , 150) ,B(3, 150) ,EMM( 150) ,XM( 150) , YM( 150) ,KOD( 150)  , 

2  IBRY( 150) ,  BV  ( 150) , MAT ( 150) ,  TMF  ( 150) ,  G  (  2  ) , CG 1 , NECH , 

3  NUMNP , NUMEL , NSEC , NE2M , MB AND ,M,MM,NN, PRIG, RIG, TMM,AR 
C 

CALL  DATA 
C 

CALL  ASSEMB 
C 

KRUN=0 
NECH=5 
DMAX=0 . 

10  CALL  SOLVE 
C 

IF (RIG- LE . PRIG)  GO  TO  35 
I F ( RIG . LE . DMAX )  GO  TO  36 
DMAX=RIG 
DO  20  1=1 , NUMEL 
20  MF IN ( I ) =MAT ( I ) 

GO  TO  36 

35  I F ( NECH . LE . 1  )  GO  TO  40 
NECH=NECH- 1 

36  PR IG=R IG 
KRUN=KRUN+ 1 

WR I TE ( 6 , 2 1 00 )  KRUN 

2100  FORMAT (  /  , '  NUMBER  OF  RUN  =',14) 

CALL  STRESS 
GO  TO  10 
C 

40  WR I TE ( 6 , 2 1 02  )  ( M , MF I N ( M ) , M= 1 , NUMEL ) 

2102  FORMAT ( 1016) 

WRITE (6,2103) DMAX 

2103  FORMAT ( '  MAX.  TORSIONAL  RIGIDITY  =' , F 1 0 . 5 ) 

STOP 

END 


SUBROUTINE  DATA 
C  THIS  SUBROUTINE  READS  AND  PRINTS  MATERIAL,  NODAL  AND 
C  ELEMENT  DATA.  IT  GENERATES  DATA  FOR  INTERMEDIATE 

C  ELEMENTS  AND  CALCULATES  THE  BAND  WIDTH  AND  NUMBER  OF 
C  EQUATIONS. 

Q  ^  ^  ^  4*  *4^  jjs  Jfc  3fC  5|c  5fC  ^  5fC  5fs  3)C  ^  4*  *4^  ^  ^  ^  jfc  Jfi  5fc  5|C  5fC  ^  *4*  '4' 

c 

COMMON  TFF( 150,30) , EFF ( 3 , 3 , 150) ,EMF(3, 150) ,NP(3, 150)  , 

1  A ( 3 . 150) ,B(3, 150) ,EMM( 150) ,XM( 150) ,  YM( 150) ,KOD( 150) , 

2  I BR Y ( 150) ,BV( 150) . MAT ( 150) ,TMF( 150) ,G(2) , CG 1 , NECH , 

3  NUMNP  ,  NUMEL , NSEC , NE2M , MBAND , M , MM , NN , PR IG , R IG , TMM , AR 
DIMENSION  HED ( 20 ) , MF IN ( 150) 
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DIMENSION  X(  150) ,Y( 150) , ISYM( 150) 

REAL  DN,DX,DY 

BASIC  DATA 

READ ( 5 , 1000)  HED,NSEC,NUMEL,NUMNP,NE2M 
1000  FORMAT)  20A4/  516) 

WRITE (6, 2000)  HED , NSEC , NUMEL , NUMNP , NE2M 
2000  FORMAT) 1H1 , 1  OX , 20A4 , //// 

1  1 H  ,  'NUMBER  OF  SECTIONS  =  '  ,  16/ 

2  1H  ,  'NUMBER  OF  ELEMENTS  =' ,  16/ 

3  1 H  ,  'NUMBER  OF  NODAL  POINTS  ='  ,  16/ 

4  1 H  ,  'NO.  OF  ELEMENTS  OF  STIFFER  MATERIAL  =' 
WRI TE ( 6 , 2005 ) 

2005  FORMAT)///, 1H  , 1  OX , '  MATERIAL  PROPERTIES  'll 
1  IX, 'MAT.  NO.'  ,5X,  '  SHEAR  MODULUS  G'  ) 

DO  10  1=1,2 
READ ( 5 , 1010)  G ( I  ) 

1010  FORMAT ( 6 X , F 1 2 . 0 ) 

10  WRITE ( 6 , 20 1  0  )  I , G 1 1  ) 

2010  FORMAT  ( 1H  , 15 , 1 1 X , F 1 0 . 3 ) 

NODAL  DATA 


L=  1 

READ) 5, 1020)  N , X ( N ) , Y ( N ) , ISYM ( N ) 

GO  TO  50 

20  READ (5,1 020  )  N , X ( N ) , Y ( N ) , ISYM ( N ) 

1020  FORMAT) I6.2F12.3, 16) 

DN=N-L 

DX= ( X ( N ) -X ( L ) )/DN 
DY=(Y(N)-Y(L) )/DN 
30  L  =  L+  1 

IF(N-L)  60,50,40 
40  X(L)=X(L-1 ) +DX 
Y ( L ) = Y ( L -  1 ) +DY 
GO  TO  30 

50  IF(NUMNP-N)  60,70.20 
60  WRITE ( 6 , 2025  )  N 

2025  FORMAT) 1  HO , '  ERROR  IN  NODAL  DATA,  NODE  ='  ,  14 ) 
CALL  EXIT 
70  WRITE ( 6 , 20 14 ) 

WRITE ( 6 , 20 1 5 ) 

WRITE (6, 2020)  (N,X(N),Y(N) ,N=1 , NUMNP) 

2014  FORMAT ( '  1 '  , 5X ,  'COMPLETE  NODAL  DATA') 

2015  FORMAT)///, 1 0X , '  NODAL  POINT  OUTPUT',/// 

1  1 H  ,'  NODE  X  COORD  Y  COORD',//) 

2020  FORMAT) 14,  2F12.2) 

ELEMENT  DATA 


ML  =  0 

80  I F ( ML . GE . NUMEL  )  GO  TO  95 

READ) 5, 1035)  M , NP ( 1 , M ) , NP ( 2 , M ) , NP ( 3 , M ) , IBRY ( M ) 
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1035  FORMAT ( 516 ) 

MM=ML+ 1 

IF(MM.GE.M)  GO  TO  90 
82  ML  1 =ML+  1 

I F ( ML  1  .  GE  .  M )  GO  TO  90 

ML2=ML+2 

MM  1 =ML -  1 

I F ( MM  1 . LE . 0 )  GO  TO  160 
DO  86  1=1,3 
NP ( I , ML2 ) =NP ( I , ML ) + 1 
NP( I , ML  1 )=NP( I ,MM1 )  +  1 
IBRY ( ML2 )  =  I B  R  Y ( ML ) 

IBRY (ML  1 )  =  I B  R  Y ( MM  1 ) 

86  CONTINUE 
ML=ML2 
GO  TO  82 
90  ML=M 
GO  TO  80 

MATERIAL  DISTRIBUTION 

95  DO  100  K= 1 , NUMEL 
100  MAT ( K )  =  1 

READ (5,1 040 )  K1.KEL1 
I F ( K 1 .GE.NE2M)  GO  TO  130 
110  READ ( 5 , 1040)  K2.KEL2 
1040  FORMAT (216) 

120  MAT ( KEL  1  )  =2 
K 1 =K 1  +  1 
KEL 1 =KEL 1+ 1 

IFIK1.LT.K2)  GO  TO  120 
KEL 1 =KEL2 

IFIK2.LT. NE2M)  GO  TO  110 
MAT ( KEL2 ) =2 
GO  TO  140 
130  MAT ( KEL  1  )  =2 
140  WRITE ( 6 , 2030 ) 

WRITE ( 6 , 2032 ) 

WRITE (6, 2035)  ( M , ( NP ( J , M ) , il  = 1 , 3 ) , MAT ( M ) , I BRY ( M  )  , 

1  M  =  1  , NUMEL) 

2030  FORMAT ( '  1 '  , 5X , '  COMPLETE  ELEMENT  DATA') 

2032  FORMAT!///,'  ELEM.  I  d  K  MAT.  BDRY' ) 

2035  FORMAT (616) 

BAND  WIDTH  AND  NUMBER  OF  EQUATIONS 
L  =  0 

DO  150  M= 1 , NUMEL 
DO  150  1=1,2 
11=1+1 

DO  150  <J  =  1 1  ,  3 
K  =  I ABS ( NP ( I ,M) -NP(d,M)  ) 

IF(K.GT.L)  L=K 
150  CONTINUE 
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MBAND=L+ 1 
NEQ=NUMNP 

WR I TE ( 6 , 2040 )  MBAND , NEQ 
I F ( MB AND . LE . 40 . AND .NEQ.LE. 150)  GO  TO  170 
WRITE ( 6 , 2050 ) 

2040  FORMAT ( ///  , 1  OX ,  '  BAND  WIDTH  =#  , 16/ 

1  10X,  ' NUMBER  OF  EQUATIONS  =',16) 

2050  FORMAT ( ///  , 1  OX , 7  PROBLEM  EXCEEDS  SPECIFIED  LIMITS') 
CALL  EXIT 
160  WR I TE ( 6 , 2060 ) 

2060  FORMAT (///,'  MM  1  IS  LESS  THAN  OR  EQUAL  TO  ZERO') 
CALL  EXIT 

165  WR I TE ( 6 , 2095 )  M 

2095  FORMAT  (  1 H  /NEGATIVE  AREA  ELEMENTS  ,  16  ) 

CALL  EXIT 

170  I F ( NSEC . LE . 1 )GO  TO  190 
MM  =  0 
NR  =  0 


NODES  ON  THE  BOUNDARY  OF  SECTIONS 

DO  180  1=1 , NUMNP 
IF ( ISYM( I ) .EQ. 1 )  GO  TO  180 
NR=NR+ 1 
KOD ( NR ) = I 
180  CONTINUE 
NN  =  NR 
RETURN 

190  MM=MBAND 
NN=NEQ- 1 
RETURN 
END 


SUBROUTINE  ASSEMB 

THIS  SUBROUTINE  TAKES  EACH  ELEMENT  IN  TURN  AND  FORMS 
THE  ELEMENT  STIFFNESS  MATRIX  AND  ASSEMBLES  THE  GLOBAL 
STIFFNESS  MATRIX. 

COMMON  TFF( 150,30) , EFF ( 3 , 3 , 150) ,EMF(3, 150) , N P ( 3 , 150) , 

1  A ( 3 , 150) ,  B(3, 150) ,EMM( 150) ,XM( 150) , YM( 150) , KOD ( 150) , 

2  IBRY ( 150) ,BV( 150) , MAT ( 150) ,  TMF  ( 150) ,G(2) ,CG1 , NECH , 

3  NUMNP, NUMEL, NSEC, NE2M, MBAND ,M, MM ,NN, PRIG, RIG, TMM, AR 

INITIALIZE 

DO  10  1=1 , NUMNP 
TMF ( I ) =0 . 0 
DO  10  J=1 , MBAND 
TFF  (  I  ,  J ) =0 . 0 
10  CONTINUE 
TMM=0 . 0 


C 
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C  ELEMENT  STIFFNESS  MATRICES 
C 

AR=(A(3, 1 )*B(2, 1 ) -A ( 2 , 1 )*B(3, 1 )  1/2.0 
CG=G(2)/G( 1 ) 

CGI =CG- 1 
GV=G (  1  ) 

C 1 =GV/2 . 

C2=GV/(4.*AR) 

DO  70  M= 1 , NUMEL 
I  =  N  P  (  1  ,  M ) 
d=NP ( 2 , M ) 

K=NP ( 3 , M ) 

A (  1  ,M)=X(K)-X(d) 

A ( 2 , M ) =X ( I )  -X ( K ) 

A ( 3 , M ) =X ( J ) -X ( I ) 

B ( 1 ,M)=Y( J)-Y(K) 

B ( 2 , M ) =Y ( K ) - Y ( I ) 

B ( 3 , M ) =Y ( I )  - Y ( d ) 

XO= ( X ( I )+X(d)+X(K) ) / 3 . 0 
YO= ( Y { I )+Y(d)+Y(K)  )/3.0 
XM ( M ) =XO 
YM(M)=YO 
A I  =  A ( 1 ,M) 

B I  =  B ( 1 ,M) 

IF ( IBRY (M) . GT . 0 )  GO  TO  30 

( 1 )  INNER  ELEMENTS 

Ad= A ( 2 , M ) 

AK= A ( 3 , M ) 

Bd  =  B  ( 2 , M ) 

BK=B ( 3 , M ) 

EMM(M)=GV*AR*(XO*XO+YO*YO) 

EMF ( 1 ,M)=C1*(AI*XO-BI*YO) 
EMF(2,M)=C1*(Ad*XO-Bd*YO) 
EMF(3,M)=C1*(AK*XO-BK*YO) 

EFF ( 1 , 1 ,M)=C2*(A1*AI+BI*BI  ) 

EFF ( 1 ,2,M)=C2*(AI*Ad+BI*Bd) 

EFF ( 1 ,3,M)=C2*(AI*AK+BI*BK) 

EFF ( 2 , 1 ,M)=EFF( 1 ,2,M) 
EFF(2,2,M)=C2*(Ad*AJ+Bd*Bd) 
EFF(2,3,M)=C2*(Ad*AK+Bd*BK) 

EFF ( 3 , 1 ,M)=EFF( 1 ,3,M) 
EFF(3,2,M)=EFF(2,3,M) 
EFF(3,3,M)=C2*(AK*AK+BK*BK) 

GO  TO  40 

(2)  BOUNDARY  ELEMENTS 

30  SL23=AI*AI+BI*BI 
EK=GV*AR/SL23 
EM=YO*AI+XO*BI 
EMM( M ) =EK*EM**2 
EMF ( 2 , M ) =EM*EK 


. 
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EMF ( 3 , M ) = -EM*EK 
EFF (2,2,M)=EK 
EFF ( 2 , 3 , M ) =-EK 
EFF(3,2,M)=-EK 
EFF ( 3 , 3 , M ) =EK 

GLOBAL  STIFFNESS  MATRIX 

40  I F  C  MAT ( M ) . GE . 2  )  GO  TO  47 
TMM=TMM+EMM(M) 

K I  =  1 

I F ( IBRY(M) . GT  .  0  )  K I =2 
DO  46  I =K I ,3 
1 1 =NP ( I , M ) 

TMF(II)=EMF(I,M)+TMF(II) 

DO  44  J  =  K 1 , 3 
dd  =  NP ( d , M ) - 1 1  + 1 
IF(dd.LE.O)  GO  TO  44 
TFF(II,dJ)=TFF(II,dd)+EFF(I,d,M) 

44  CONTINUE 

46  CONTINUE 
GO  TO  70 

47  TMM=TMM+CG*EMM ( M ) 

K I  =  1 

I F  < IBRY(M) . GT . 0  )  K I =2 
DO  49  I=KI ,3 
1 1 =NP ( I , M ) 

TMF (II) =CG*EMF (I,M)+TMF(II) 

DO  48  d  =  KI  ,3 

dd=NP(d,M)-II+1 

IF(dd.LE.O)  GO  TO  48 

TFF  ( 1 1 , dd ) =TFF (II, dd ) +CG*EFF ( I , d,M) 

48  CONTINUE 

49  CONTINUE 
70  CONTINUE 

RETURN 

2150  FORMAT ( 10F7 . 3) 

END 

SUBROUTINE  SOLVE 

THIS  SUBROUTINE  SOLVES  THE  SET  OF  EQUATIONS  FOR 
THE  VALUES  OF  WARPING  FUNCTION  AT  THE  NODAL  POINTS 
AND  CALCULATES  THE  TORSIONAL  RIGIDITY. 

COMMON  TFF ( 150,30) , EFF (3, 3, 150) , EMF (3, 150) ,NP( 3, 150) , 

1  A ( 3  ,  150) ,B(3, 150) , EMM ( 150) ,XM< 150) ,YM( 150) , KOD ( 150) , 

2  I BRY  (  150) ,BV(  150) , MAT ( 150) , TMF ( 150) ,G(2) , CGI  ,  NECH  , 

3  NUMNP , NUMEL, NSEC, NE2M,MBAND,M, MM, NN, PRIG, RIG, TMM.AR 
DIMENSION  AMI  150,30) 

IFINSEC.LE.  1  )  GO  TO  16 
IF(MM.GT.O)  GO  TO  10 
DO  3  1=1 , NN 
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DO  3  d=1 , MBAND 

3  AM ( I , J ) =0 . 0 
1=0 

4  1=1+1 

IF ( I .GT.NN)  GO  TO  14 
1 1 =KOD ( I ) 

BV ( I ) =-TMF (II) 

KOUNT  =  0 
J  =  I  -  1 

5  d  =  d+ 1 

J J  =  KOD ( J ) - 1 1  + 1 

IF(dd.GT. MBAND)  GO  TO  4 

L= J- 1+ 1 

KOUNT  =  KOUNT+ 1 

AM ( I , L ) =TFF ( 1 1 , J J ) 

IF ( MM. LT. KOUNT)  MM=KOUNT 
IF(J.LT.NN)  GO  TO  5 

7  1=1+1 

I F ( I . GT.NN)  GO  TO  14 
1 1 =KOD (  I  ) 

BV ( I ) = - TMF (II) 

DO  8  d=I,NN 
J J  =  KOD ( J ) - 1 1  + 1 
IF(JJ.GT. MBAND)  GO  TO  9 
L= J- 1+1 

8  AM ( I , L ) =TFF ( 1 1 , J J ) 

9  GO  TO  7 

10  NL=NN-MM+ 1 
DO  12  1=1 , NN 
1 1 =KOD ( I ) 

BV ( I ) = - TMF (II)  • 

DO  11  d  =  I ,  NN 
L  =  J  - 1  + 1 

IF(L.GT.MM)  GO  TO  12 
J J  =  KOD ( J ) - 1 1  +  1 
IF(dd.GT. MBAND)  GO  TO  12 

11  AM ( I , L ) =TFF ( 1 1 , dd  ) 

12  CONTINUE 
14  NL=NN-MM+ 1 

NM  =  NN-  1 

IF(NM)  110,90,30 
16  NL=NN-MM+ 1 
NM  =  NN-  1 

STORE  TFF  AND  TMF  IN  AM  AND  BV  AND  APPLY  DISPLACEMENT 
BOUNDARY  CONDITION. 

DO  18  K  = 1 , NN 
K 1 =K+ 1 

18  BV ( K ) = - TMF ( K  1  ) 

DO  20  d= 1 , MM 
DO  20  1=1 , NN 
11=1+1 

20  AM ( I , d ) =TFF ( 1 1 , d ) 


ooo  ooo  ooo 


72 


TRIANGULARIZE  AM  AND  REDUCE  BV 


30  MR=MM 

DO  52  N= 1 , NM 

I F ( AM ( N , 1 ) .GT. 0.0001 )  GO  TO  35 
AM ( N ,  1  )  =  1  . 

BV ( N ) =0 . 

DO  32  KC=2 , MM 
32  AM ( N , KC ) =0 . 

GO  TO  52 
35  BN=BV ( N ) 

BV ( N ) =BN/AM ( N , 1 ) 

IF(N.GT.NL)  MR=NN-N+ 1 
DO  50  L=2 , MR 

I F ( AM ( N , L ) .EQ.0.0)  GO  TO  50 
C  =  AM ( N , L ) /AM ( N ,  1  ) 

I =N+L -  1 
J  =  0 

DO  40  K=L , MR 
d  =  d+  1 

40  AM ( I , J ) = AM ( I , d ) -C*AM ( N , K ) 

B  V ( I ) =  B  V ( I ) -C*BN 
AM ( N , L  )  =C 
50  CONTINUE 
52  CONTINUE 

BACK  SUBSTITUTE 

I  =NN 

55  BV(NN)=BV(NN)/AM(NN, 1 ) 

DO  60  N= 1 , NM 
1  =  1-1 

IF(N.LT.MM)  MR=N+ 1 
DO  60  d=2 , MR 
K=  1  + J-  1 

60  BV ( I ) =  B V ( I ) - AM ( I ,  d ) *BV ( K ) 

* 

TORSIONAL  RIGIDITY 
T  =  0 . 0 

I F ( NSEC . GT . 1  )  GO  TO  82 
70  DO  80  11=1 ,NN 
I d  =  1 1  + 1 

T  =T+TMF (Id) *BV (II) 

80  CONTINUE 
GO  TO  88 

82  DO  85  11=1 , NN 
I d=KOD (II) 

85  T  =T+TMF (Id) *BV (II) 

88  R I G  = ( TMM+T ) *NSEC 
WRITE ( 6 , 2070  )  RIG 
RETURN 

90  B V ( 1 ) =  B V ( 1 ) /AM ( 1,1) 
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dP=K0D (  1  ) 

RIG=(TMM+TMF(dP)*BV( 1 ) )*NSEC 
WRITE ( 6 , 2070  )  RIG 

2070  FORMAT (///,'  TORSIONAL  RIGIDITY  ='  , F 1 2 . 6 ) 

RETURN 

110  WRITE ( 6 , 30 1 0 ) 

3010  FORMAT!' ELEMENTS  ARE  NOT  DIVIDED  PROPERLY') 

CALL  EXIT 

3000  FORMAT! 1H  ,'ZERO  OR  NEGATIVE  ELEMENT  ON  MAIN  DIAGONAL 
1  OF  TRIANGULARIZED  MATRIX' ,  15) 

END 

C 

SUBROUTINE  STRESS 

THIS  SUBROUTINE  COMPUTES  STRESSES  IN  EVERY  ELEMENTS, 
COMPARE  THEIR  MAGNITUDES,  SELECT  PAIRS  OF  ELEMENTS 
FOR  INTERCHANGE  AND  MODIFY  THE  STIFFNESS  MATRIX. 

COMMON  TFF!  150,30) , EF  F ( 3 , 3 , 150) ,EMF(3, 150) , NP ( 3  ,  150)  , 

1  A ( 3 , 150) ,B(3, 150) ,EMM( 150) ,XM( 150) ,YM( 150) , KOD ( 150) , 

2  IBRY!  150) ,BV(  150) .MAT! 150) , TMF ( 1 50 ) , G ( 2 ) , CGI , NECH , 

3  NUMNP,NUMEL,NSEC,NE2M,MBAND,M,MM,NN, PRIG,  RIG, TMM,AR 
DIMENSION  W( 150) , SQT! 150) 

DO  10  1=1, NUMNP 

10  W( I ) =0 . 

I F ( NSEC . LE . 1 )  GO  TO  30 
DO  20  d=1 ,NN 
KP=KOD(d) 

20  W ( KP ) =BV ( d ) 

GO  TO  40 

SHEAR  STRESSES 

30  DO  35  J= 1 , NN 
d1=d+1 

35  W ( 0 1 )=BV(d) 

40  DO  60  K= 1 , NUMEL 
I  1 =NP ( 1 , K ) 

1 2  =  NP ( 2 , K ) 

I3=NP ( 3 , K  ) 

IF! IBRY(K) . NE . 0  )  GO  TO  50 

TXZ=( -YM(K)  +  (B( 1 ,K)*W( 1 1 )+B(2,K)*W( 12 ) +B ( 3 , K ) *W( 13) )/ 
1  ( 2 . *AR ) ) 

TYZ=(XM(K)  +  ( A! 1 ,K)*W( 1 1 )+A(2,K)*W( 12 )+A(  3 ,K )*W(  13)  )/ 

1  ( 2 . *AR )  ) 

GO  TO  55 
50  AK=A ( 1 ,K) 

BK=B ( 1 ,K) 

SL23=AK*AK+BK*BK 

TXZ=(AK/SL23)*(-AK*YM(K)-BK*XM(K)-W(I2)+W(I3) ) 
TYZ=(BK/SL23)*(AK*YM(K)-BK*XM(K)+W(I2)-W(I3) ) 

55  SQT(K)=TXZ*TXZ+TYZ*TYZ 
60  CONTINUE 
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SELECT  PAIRS  OF  ELEMENTS 
KAL=  1 

65  TMIN=SQT ( 1 ) 

TMAX=TMIN 
NDN=  1 
NUP=  1 

70  DO  80  M=2 , NUMEL 

IF (MAT(M) . GE . 2  )  GO  TO  75 
IF ( SQT (M) .LE.TMAX)  GO  TO  80 
TMAX=SQT ( M ) 

NUP  =  M 
GO  TO  80 

75  I F ( SQT ( M ) .GE.TMIN)  GO  TO  80 
TMIN=$QT ( M ) 

NDN  =  M 

80  CONTINUE 

IF(TMIN.GE.TMAX)  GO  TO  100 
K I  =  1 

MAT ( NUP ) =2 
MAT ( NDN )  =  1 

IF ( IBRY(NUP) .NE.O)  KI=2 

MODIFY  THE  SET  OF  EQUATIONS 

DO  95  I =  K I  ,3 
1 1 =NP ( I , NUP ) 

TMF ( II )=TMF ( II ) +CG 1 *EMF ( I , NUP ) 

DO  93  d  =  KI  ,3 
J J  =  NP ( J , NUP ) - 1 1  + 1 
IF(dd.LE.O)  GO  TO  93 
TFF(II , J J ) =TFF ( 1 1 , dd)+CG1*EFFd f  J , NUP ) 

93  CONTINUE 
95  CONTINUE 

K I  =  1 

IF(  IBRY(NDN) .NE.O)  KI  =  2 
DO  99  I  =  K I ,3 
1 1  =  NP ( I , NDN ) 

TMF ( 1 1 ) =  TMF ( 1 1 )-CG1*EMF( I , NDN ) 

DO  94  J=KI ,3 
J J  =  NP ( J , NDN ) -  1 1  +  1 
IF(dd.LE.O)  GO  TO  94 

TFF(II , dd ) =TFF ( 1 1 ,dd)-CG1*EFF(I , d , NDN ) 

94  CONTINUE 
99  CONTINUE 

TMM  =  TMM+CG 1  * ( EMM ( NUP ) -  EMM ( NDN ) ) 

WRITE ( 6 , 2805 )  NDN, NUP 
IF(KAL.GE.NECH)  GO  TO  100 
KAL=KAL+ 1 

2805  FORMATS  ELEMENT  CHANGED  #'  ,16/  AND  #',I6) 
GO  TO  65 
100  RETURN 
END 


